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FIRST YEAR PROGRESS REPORT 


This report covers technical progress during the first year of the NASA Space Physics 
Theory contract “The Structure and Dynamics of the Solar Corona,” NAS5-96081, between 
NASA and Science Applications International Corporation, and covers the period July 26, 1996 to 
July 16, 1997. Under this contract SAIC, the University of California, Irvine (UCI), and the Jet 
Propulsion Laboratory (JPL), have conducted research into theoretical modeling of active regions, 
the solar corona, and the inner heliosphere, using the MHD model. During the period covered by 
this report we have published 26 articles in the scientific literature. These publications are listed in 
Section 4 of this report. In the Appendix we have attached reprints of selected articles. 

1. INTRODUCTION 

An important goal of NASA’s space physics program is to develop a quantitative 
understanding of the Earth’s plasma environment. The Sun, as the driver of terrestrial 
geomagnetic disturbances, has an all-important influence on the Earth. This Sun-Earth connection, 
dubbed “space weather,” has recently been recognized as an important activity for the space science 
community. 

Our program centers around the theoretical modeling of active regions, the solar corona, and 
the inner heliosphere, using the MHD model. This capability allows us to compare our model 
quantitatively with many observations, including interplanetary scintillation (IPS) and in situ 
spacecraft measurements of the solar wind, and coronal emission (e.g., white light, radio, SOHO 
EIT and UVCS). 


2. ACHIEVEMENTS 

In this section we summarize the accomplishments made by our group during the first year of 
our Space Physics Theory Program contract. The descriptions are primarily intended to illustrate 
our principal results. A full account can be found in the referenced publications. 

2.1. Modeling the Large-Scale Structure of the Corona and Inner Heliosphere 

Our description of the solar corona and inner heliosphere is based on numerical solutions of 
the MHD equations. We have extended the state-of-the-art in coronal modeling in three principal 
areas: we have brought in geometrical realism (by extending ID and 2D models to 2D and 3D 
models, respectively); we have used more sophisticated physics, including an improved 
description of energy flow; and we have used solar observations, principally the measured 
photospheric magnetic field, to make it possible to directly compare our calculations with solar 
observations. 

2JJ. An Improved MHD Model of the Solar Corona and Inner Heliosphere 

Our modeling of the global properties of the solar corona relies on the MHD model to 
describe the interaction of the solar wind with coronal magnetic fields. In order to simplify the 
description, we initially used a “polytropic model,” in which an adiabatic energy equation with a 
reduced poly tropic index y (i.e., smaller than 5/3) is used (Parker 1963). This is a crude way of 
modeling the complicated thermodynamics in the corona with a simple energy equation. The 
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primary motivation for using a reduced y is the fact that the temperature in the corona does not 
vary substantially, since the limit y — » 1 corresponds to an isothermal plasma. A typical choice 
is y ~ 105 - Comparisons of our results with coronal observations indicate that while this model 
matches many features of the corona, it is not accurate enough to quantitatively reproduce the 
properties of the corona and solar wind. In particular, this simple model fails to reproduce the fast 
(~ 800 km/s) and slow (~ 400 km/s) solar wind streams that are measured at 1 AU, nor does it 

reproduce the contrast in density and temperature that is observed between streamers and coronal 
holes. 

We have improved this aspect of the formulation by modeling in detail the physical 
mechanisms that describe the transport of energy in the corona and solar wind. One-dimensional 
MHD models have been quite successful, despite their obvious geometrical limitations, in 
describing this interaction and in comparisons with spacecraft solar wind measurements (Withbroe 
1988; Habbal et al. 1995). We have improved the energy equation in our model to include the 
effects of parallel thermal conduction, radiation loss, parameterized coronal heating, and Alfven 
wave acceleration. Our improved model is based on the following MHD equations: 

Air 
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S = ~ V q - n e n p Q(D + H ch + H d + D , (7) 

where B is the magnetic field, J is the current density, E is the electric field, p, v, p, and T are the 
plasma mass density, velocity, pressure, and temperature, and the wave pressure p w represents the 

acceleration due to Alfven waves. The gravitational acceleration is g, y = 5/3 is the ratio of 
specific heats, rj is the resistivity, v is the viscosity, // c h is the coronal heating source, D is the 
Alfven wave dissipation term, n e and n p are the electron and proton density, and Q{T) is the 
radiation loss function (Rosner et al 1978). The term H d = pj 2 + vVv:Vv represents heating 
due to viscous and resistive dissipation. In the collisional regime (below ~ 10 R s ), the heat flux is 
q = -/c N bb-Vr, where b is the unit vector along B, and Kj, = 9 x 10 S * 7 T 5/2 is the Spitzer value 
of the parallel thermal conductivity. In the collisionless regime (beyond - 10/y, the heat flux is 
given by q = an e kT\ , where a is a parameter (Hollweg 1978). Since it is presently not known 
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in detail what heats the solar corona, the coronal heating source H c h is a parameterized function. A 
typical form is 

H ch = # o (0)exp [-(r-/y/;i(0)] , (8) 

where H Q (6) expresses the latitudinal variation of the volumetric heating, and A(0) expresses the 
latitudinal variation of the heating function scale length. [In practice, the variation is expressed in 
terms of the magnetic topology (i.e., a proxy for the open and closed field regions) rather than the 
latitude 0.] Note that the simplified polytropic model is obtained by setting S = 0 in Eq. (6), 
P w - 0 in Eq. (5), and y = 1-05. 

Since the acceleration of the solar wind by Alfven waves occurs on a spatial and time scale 
that is below the spatial and time resolution of our global numerical model, the wave pressure p w is 
evolved using an equation for the time-space averaged Alfven wave energy density e 
(Jacques 1977), 

^+V-F = y-V Pw -D , (9) 

where F = (- v + \ A )e is the Alfven wave energy flux, v A = BI^Anp is the Alfven speed, and 
Pw - 2 £ ’ The Alfven wave velocity is = ±1)^; in a multi-dimensional implementation, it is 
necessary to transport two Alfven wave fields (waves parallel and antiparallel to B), which are 
combined to give e. The Alfven wave energy density e is related to the space-time average of the 
fluctuating component of the magnetic field SB by e = <8B 2 >!Ak. The dissipation term D 
expresses the nonlinear dissipation of Alfven waves in interplanetary space and is modeled 
phenomenologically (Hollweg 1978). 

We use the following boundary conditions (Mikic & Linker 1994; Linker et al. 1996; Linker 
& Mikic 1997). The radial magnetic field is specified at the solar surface r = R s (e.g., from 
synoptic magnetic field observations, or from full-disk magnetograms). The boundary conditions 
on the velocity are determined from the characteristic equations along B. The boundary r = R s is 
chosen to be at the top of the transition region, at a given temperature (say T 0 = 500,000°K). The 
density at r = R s is determined by balancing radiation loss, thermal conduction, and heating 
within the chromosphere and transition region (Withbroe 1988). In this formulation, the only 
boundary conditions required from observations at the base of the corona r = R s are on the radial 
magnetic field. [In the polytropic model, in contrast, we specify a uniform density in the open- 
field regions and a uniform temperature at the base of the corona.] 

We have developed a three-dimensional code to solve MHD equations (l)-(9) in spherical 
coordinates (r,0,0). This code is described in part by Mikic and Linker (1994) and by Lionello, 
Mikic, and Schnack (1998). This code has been used extensively to model the 2D and 3D corona, 
including the solar wind, differential solar rotation, simulation of the interplanetary medium (from 
the Sun to 1 AU), the effect of emerging flux, coronal mass ejections, and the long-term evolution 
of the solar corona and heliospheric current sheet. In the following sections we describe the 
application of this code to several problems of interest in large-scale coronal physics. 

This improved model makes it possible to study the solar wind from its origin in the low 
corona to its expansion into interplanetary space. It is presently being used in 2D, and will soon be 
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implemented in 3D. We still routinely employ the polytropic model (currently the only model we 
use in 3D) when a qualitative description of the corona is sufficient, since the full model requires 
considerably increased computational resources. With the 3D improved model it will be possible 
to directly compare model results with in situ spacecraft measurements. 

Even with these improvements to the energy equation, it must be recognized that a single- 
fluid description (inherent in the MHD model) is still a considerable approximation to the state of 
the corona. In particular, recent SOHO observations imply that the electron temperature is 
considerably lower in the corona than the ion temperature. One-dimensional models (e.g., Li, 
Esser, & Habbal 1997) have extended the theory to multiple fluids. Our present philosophy is to 
include the effects of multi-dimensional geometry in a single-fluid model first, since this effect has 
not been explored fully. Eventually our formalism can be extended to include multiple fluids. 

2J.2. Comparison with Eclipse Observations 

We have by now established a “tradition” of using our 3D MHD model to predict the state of 
the corona during forthcoming total solar eclipses. Our first attempt was a comparison that was 
performed subsequent to the eclipse of November 3, 1994. (Our eclipse comparisons can be seen 
on our Web page, http://iris023.saic.com:80(X)/corona/modeling.html.) Since then, we have made 
two predictions , before the actual eclipse date, using magnetic field data from the previous solar 
rotation. Our success has been reassuring, primarily because these eclipses occurred close to solar 
minimum, when the large-scale structure of the Sun changes slowly between solar rotations. In 
Figure 1 we compare our comparisons and predictions to observations of the eclipses of November 
3, 1994, observed in Chile, October 24, 1995, observed in Vietnam, and March 9, 1997, observed 
in Siberia. 

To perform these calculations, we used Kitt Peak National Observatory and Wilcox Solar 
Observatory synoptic magnetic field maps to specify the radial magnetic field at the photosphere as 
a boundary condition for our 3D MHD solutions to compute the state of the coronal plasma, 
including the density and magnetic field (Mikic & Linker 1996). The simulated polarization 
brightness from our simulation is computed by integrating the electron density along the line of 
sight in the plane of the sky. Comparisons with Mauna Loa MK3 coronagraph observations 
during the solar rotation surrounding the eclipse have confirmed that the basic large-scale three- 
dimensional structure of the streamer belt has been captured in the model. 

2. 1.3. M HD Modeling for the Uhsses Fast Latitude Scan 

We have also applied our model to Ulysses observations during the fast latitude scan 
(Feb.~Apr., 1995). We used the MHD model to deduce the solar origin of plasma observed at 
Ulysses. HCS crossings from the MHD model were compared with those from the source-surface 
model and Ulysses measurements. Figure 2 shows that the polarity of the magnetic field predicted 
by the MHD model generally matches Ulysses observations. It is apparent that the Ulysses data 
shows a general trend: the fast wind comes from deeper within coronal holes than the slow wind . 
Full details are given by Neugebauer et al (1998). 
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(MHD Model) 


Polarization Brightness 
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October 24, 1995 
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Figure 1. Comparison of MHD computations of the solar corona with total solar eclipse 
observations. The 1994 eclipse image is courtesy of the the High Altitude Observatory, NCAR, 
Boulder, Colorado, USA. NCAR is sponsored by NSF. The 1995 eclipse image is courtesy of F. 
Diego and S. Koutchmy. The 1997 eclipse image is courtesy of the Eclipse Team of Meisei 
University, lead by Professor Eijiro Hiei of Meisei University and the National Astronomical 
Observatory of Japan. 
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Figure 3. Emergence of a current-carrying coronal loop from 
below the photosphere. Note that the field lines twist around the 
axis of the loop as they traverse from one footpoint to the other. 
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2.1.4. Modeling the Solar Wind 

The improvements to the energy equation described in Section 2.1.1 were implemented 
primarily to model the properties of the solar wind. We have performed extensive modeling in 2D 
(axisymmetric) geometry in order to self-consistently model the solar wind from its origins in the 
low corona to its expansion into interplanetary space. In order to match observed fast and slow 
wind velocities, mass fluxes, and densities at 1 AU, the coronal heating scale length (A in Eq. 9) 
needs to be shorter in the streamer belt than near the poles, consistent with Withbroe (1988). We 
have also found that it is necessary to specify an Alfven wave flux in order to match the fast wind 
speed. Our results show promising matches with generic in situ observations of the fast and slow 
solar wind, as well as the observed properties in the low corona (density and temperature contrasts 
between the streamer belt and coronal holes). This model is directly extensible to 3D geometry. 

2*1* 5- Correlation of Flux Tube Expansion with Solar Wind Velocity 

Past observational work has demonstrated an anti-correlation between coronal magnetic flux- 
tube expansion and solar wind speed (Wang & Sheeley 1990). We have compared solar wind 
velocities from Ulysses/SWOOPS to magnetic expansion factors computed from a 3D MHD model 
for Carrington rotation 1892 (during the Ulysses “fast latitude scan”; See Sec. 2.1.3). The 
comparison was made for a period during which Ulysses saw both fast and slow solar wind. We 
found that the MHD model showed slower wind from an equatorial coronal hole with a larger 
expansion factor than from a coronal hole with a smaller expansion factor, consistent with previous 
observations. However, in general, we found the correlation between solar wind speed and 
expansion factor to be very weak (Liewer et al. 1996). 

2.2. Coronal Mass Ejections 

2.2.1. E ruption of 3D Streamers bv Photospheric Shear 

As an extension of our study of the disruption of axisymmetric arcades (Mikic & Linker 
1994) and helmet streamers (Linker & Mikic 1995), we studied the evolution of helmet streamers 
in 3D geometry. We found that a similar disruption occurs in 3D geometry, and that the ejected 
plasmoid has a large-scale longitudinal extent, as seen in LASCO observations in which CMEs 
appear to leave the Sun on both limbs simultaneously. A discussion of the role of photospheric 
shear in initiating CMEs is given by Mikic & Linker (1997). The propagation of the ejection to 
1 AU is discussed by Linker & Mikic (1997). 

2.3. Three-Dimensional Active Region Modeling 

23.1. Emergence o f Twisted Coronal Loops 

In a previous study we investigated the properties of coronal loops that were formed 
dynamically by twisting motions in the photosphere (Van Hoven, Mok, & Mikic 1995). Using the 
new emerging flux capability in our 3D Cartesian MHD code, we have studied the emergence of 
twisted (i.e., current-carrying) magnetic flux in active regions (Leka et al. 1996). Since we do not 
wish to simulate the dynamics of the convection zone, we take a phenomenological approach to 
flux emergence: we specify time-dependent profiles of the normal magnetic field, B z and the 


7 



normal current density, J z as boundary conditions at the base of the corona (z = 0). The profiles 
B z and J z are specified by applying a tangential electric field at the boundary. 

This formalism gives us considerable flexibility in specifying the properties of the emerging 
magnetic flux and current in order to model the observed magnetic structures in the solar 
atmosphere. (The profiles B z and J z could be inferred from vector magnetograph observations.) 
By appropriately specifying B z and J z we have modeled the emergence of twisted magnetic field 
loops (Mok, Van Hoven, & Mikic 1997). Figure 3 shows the emergence of a single magnetic 
loop. Note that the field lines are twisted, showing that the loop carries current. The loop is non- 
planar and truly three-dimensional, displaying an S-shape when viewed from above. 

23,2. Comparison with Radio Observations 

In a collaboration with J. Lee (U. Maryland) we compared a force-free coronal magnetic field 
extrapolation (Mikic & McClymont 1994; Mikic, Linker, & Schnack 1996), determined from 
vector magnetograph observations of active region 6615 in May 1991, with observed radio 
emission. We found that force-free fields match the radio emission better than potential fields do 
(Lee et al. 1998a, c). The radio data and magnetic field model were used to deduce the density in 
the active region, and to compare it to observations (Lee et al. 1998b). 

2.4. Coronal Heating 

2.4.1. Coronal Heating and Magnetic Energy Release 

As part of our study of coronal heating, we have performed 3D numerical simulations of the 
boundary-driven quasi-static magnetic-dissipation model proposed by Parker (1972, 1983) at 
moderately high values of the Lundquist number S. We have concentrated on the dynamics of the 
electric current and the formation of intense localized current filaments whose dissipation via 
magnetic reconnection provides an ohmic heating source. We have shown that the coronal plasma 
response to quasi-static photospheric motions is to produce current filaments via an ideal MHD 
instability. The line-tied coalescence instability is the most promising model for the observed 
instability. The detailed results are discussed by Hendrix (1996) and Hendrix & Van Hoven 
(1997). A detailed description of the reconnection that occurs in the current sheets is described by 
Schnack, Mikic, & Hendrix (1998). 

2.4.2. Coronal Heating by Dissipation of Energy in Coronal Is>nps 

The kink instability of coronal loops plays a fundamental role in the heating of the solar 
corona. Loops can become linearly unstable, resulting in concentration of the current density in the 
nonlinear phase of the instability (Einaudi, Lionello, and Velli 1997). The dissipation of these 
intense current layers is an effective mechanism for converting magnetic energy into thermal 
energy. We have studied different classes of field models (Velli, Lionello, and Einaudi 1997) 
using an MHD code in cylindrical coordinates (Lionello, Mikic, & Schnack 1998). The nonlinear 
evolution of the equilibria shows that loops formed by the slow twisting of magnetic field lines 
(Mikic, Schnack, & Van Hoven 1990) are more likely to liberate large amounts of energy 
(Lionello, Velli, Einaudi, & Mikic 1998). The destabilizing role of line-tying in the nonlinear 
phase of the instabilities has also been addressed (Lionello, Schnack, Einaudi, & Velli 1998). 
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2.5. Massively Parallel Computing 

To improve the efficiency of our computations, and, in particular, to extend our ability to do 
high-resolution runs, we have investigated how to implement the spherical 3D MHD code on the 
high performance parallel supercomputers at JPL. The algorithm for the parallel implementation 
has been designed- A two-dimensional decomposition of the three-dimensional computation 
domain will be used. Only the (/*— G) plane will be divided among the processing elements (PEs) so 
that each PE has a range of r and £?, but all 0. The Fourier transforms in <j) can be done without 
interprocessor communication. At a meeting at SAIC between SAIC, JPL and Caltech researchers, 
it was decided that, in order to prevent a major code re-write in the future, the code will be re- 
written in a modern language (High Performance Fortran or Fortran 90) before beginning the 
parallel implementation. The code will be implemented on the new JPL parallel supercomputer, a 
Hewlett Packard SPP2000, which uses a PA-8000 chip. The building block of the SPP2000 
computer system is a hypemode. Each hypemode is a 16-CPU computer by itself; 4 hypemodes 
joined together form a supemode. The final configuration for the new machine will be either be a 
two 128-PE system, or one 256-PE system. The peak speed is 720 Mflops per PE. Sequential 
code can be run on the machine using global shared memory. At present, the maximum global 
shared memory is among 64 PEs only. In the future Caltech will help HP to make the global 
shared memory work on all 256 PEs. The machine supports both Fortran 90 and High 
Performance Fortran. In order to make our 3D MHD code an efficient parallel code, at present we 
are leaning to implementing interprocessor communication using a message passing interface 
(MPI). 
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ABSTRACT 

The mechanism of the dynamic emergence of current-carrying magnetic loops into the corona is examined by 
three-dimensional MHD simulations. By simultaneously modeling the spacetime profiles of the normal components 
of the emerging magnetic field and current density on the photosphere, we demonstrate that this process can 
qualitatively reproduce observations that show the emergence of a helically twisted magnetic structure with a 
suitable field-current combination. The flux-tube structure rises into the initially nearly field-free corona and 
gradually relaxes into a nearly force-free, steady state, magnetic loop. 

Subject headings: magnetic fields — MHD — plasmas — Sun: corona — Sun: magnetic fields 


1. INTRODUCTION 

Coronal loops in the solar atmosphere have been widely 
observed since the Sky lab era (Foukal 1976; Krieger 1978; Van 
Hoven 1981). Their interpretation as plasma loops collimated 
by active region magnetic fields was later reinforced by Yohkoh 
observations (Klimchuk et al. 1992). Because of their long 
duration in the solar atmosphere, they are believed to be in 
stable equilibrium with their environment. Consequently, nu- 
merous investigations (using a cylindrical model) have been 
devoted to studying their physical properties, such as their in- 
ternal geometry, force balance, and magnetohydrodynamic sta- 
bility (Raadu 1972; Giachetti et al. 1977; Hood & Priest 1979; 
Einaudi & Van Hoven 1983; Mikic. Schnack, & Van Hoven 
1980). 

However, there have been few studies of how these magnetic 
loops, which carry electric current through the corona, come 
into existence or of the effects of their quasi-toroidal geometry. 
In a previous investigation, we addressed some of these issues 
and reported that these loops can be formed in one way by 
photospheric plasma vortex motions driven by high-& (ratio of 
fluid to magnetic forces) solar-surface convection (Van Hoven, 
Mok, & Mikic 1995). Vortical motions at the poles twist the 
in situ bipolar field lines at the base and induce an azimuthal 
field component into the field structure as well as a toroidal 
current along the to-be-formed loop. 

However, recent observations also support a second scenario 
in which magnetic flux already carries electric current at the 
time of emergence into the corona (Leka et al. 1996). The 
physical process of emerging current-carrying flux into a nearly 
field-free corona is fundamentally different from twisting ex- 
isting field lines. While the latter does not change the total flux 
at the base, the former injects new flux through the base. How 
the plasma-field in the corona responses to this injection is 
reported in this Letter. To simulate the observation, we use a 
numerical model to describe the dynamic formation of magnetic 
loops within the framework of three-dimensional MHD. The 
results provide the first realization of the simultaneous emer- 
gence of toroidal fields and currents, in contrast to the work 
of McClymont & Mikic (1994) in which only current was 


injected. We succeed in representing one physical mechanism 
for the often observed emergence of flux into the solar corona. 

2. PHYSICAL AND COMPUTATIONAL MODELS 

We model this type of flux emergence by starting with a 
small background magnetic field embedded in the coronal 
plasma. Magnetic field and current density are then injected 
into a Cartesian computational domain at the base. It is argued 
that magnetic flux, coexisting with axial current density, is in 
the form of a flux rope created below the surface in a high-/? 
plasma. As buoyancy brings the loop to the surface, the com- 
bined, self-consistent, field/current simply appears on the pho- 
tosphere. From a computational point of view, we are treating 
this process as a phenomenological problem. That is, we are 
interested in understanding how the field/current develops in 
the corona as a response to the evolving conditions on the 
photospheric surface. The physical process under the surface 
is not fully understood (Fan, Fisher, & McClymont 1994; Lo- 
ngcope, Fisher, & Arendt 1996) and is not part of our com- 
putational domain. 

Our approach in this work is justified because the photo- 
sphere is a high-/? plasma, which can control how the flux 
should emerge, and is a boundary of the computational domain, 
which is strictly regarded as the low-/? corona. We merely 
model the magnetic field and the current density at this bound- 
ary, according to observations, to provide a driving term in the 
dynamic equations. To do this, we specify the normal com- 
ponents of the magnetic field and the current density at this 
surface to simulate those of the emerging flux rope. Further 
specification on the boundary, such as the tangential compo- 
nents of either quantity, will result in overspecification and 
inconsistency in our algorithm. 

We describe the coronal dynamics with a reduced set of the 
full MHD equations. The main nonstandard assumption used 
is that of constant density and pressure as justified (partially a 
posteriori) by Ortolani & Schnack (1993) for nearly force-free 
fields in the absence of significant energy transport. The phys- 
ical motivation behind this assumption is that the corona is a 
low-/3 plasma and thus density-variation and plasma-pressure 
effects are most likely negligible in the presence of the dom- 
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Fig. 1. — Time profiles of the emerging B n . and J, h at the center of the poles 
at the base, along with five actual peak values of B , and J, at the ape. r. 

inant magnetic forces. Under these approximations , we are left 
with a system governed by two coupled differential equations 
describing magnetic induction and fluid momentum transport: 


dA/dt — r X B — riJ, ( 1 ) 


pidvldt + v * Vv) = 7 x 5 + vpV 2 v, ( 2 ) 

where A is the vector potential, v is the plasma flow velocity 
(maintained below 3% of the Alfven speed throughout the ev- 
olution), p is the plasma density, 17 and v are the resistivity and 
viscosity, respectively, written in the customary dimensionless 
units of MHD (Ortolani & Schnack 1993). The magnetic field 
B and current density / are given by 

B~VxA and J-VxB. (3) 


The implementation of the emerging magnetic field and cur- 
rent density begins with a pair of predefined spacetime profiles 
for B 0 . and 7 0; , specified only at the base. Since the vector 
potential A is the actual dynamic variable, with B and J being 
postcomputed quantities, the desired values of B. and 7. at the 
boundary, i.e., fi 0 . and 7 0 _, must be integrated to obtain the 
corresponding evolving value of A. Because of the staggered 
grids used for A in the finite-differencing scheme, only its 
tangential components A, are imposed as a time-varying bound- 
ary condition. To advance A, at the base, we first evaluate the 
tangential electric field E , contributed by B 0 . and7 0 ... In general, 
E, can be expressed as V,x zip — where the first part is 
induced by the emerging B 0 . and the second is determined by 
7 0 .. The induced part \p can be obtained by solving the equation 

V, 2 ^ = dB 0 Jdt. (4) 


To achieve the desired current density 7 0 . at the base, an 
electric potential <t> is also applied at the boundary. The purpose 
is to create a potential difference between the two ends of each 


field line in order to drive a current. Effectively, each field line 
through the corona is a circuit with an equivalent battery con- 
necting its two ends beneath the surface. Since we are only 
modeling the intended electric current 7 0 . phenomenologically, 
we are not concerned with the actual mechanism below the 
surface that drives this current. In order to make the actual 
computed 7. approach the desired J Q . , <t> is chosen to satisfy the 
equation 

d<t>/dt = (7 0; - 7.)/r,., (5) 

where t c is a time constant, chosen to be ~ 0. 1 of the rise time 
of the emerging B., so that the current will closely follow the 
field. This source-term algorithm has also been used by 
McClymont & Mikid (1994; McClymont, Jiao, & Mikic 1997) 
in the presence of fixed B 0: sources. 

Finally, the tangential vector potential A, at the base evolves 
in time according to Maxwell’s equation 

dAJdt = -E,. (6) 

With this time-dependent boundary condition for A, at the base, 
equation ( 1 ) is advanced with additional conditions on the other 
five surfaces where B n and E, vanish. The momentum equation 
(2) is advanced with the base boundary condition that Ohm’s 
law is satisfied, 

v = E x B/B\ (7) 

and v = 0 on the other five surfaces. 

The magnetic loop is taken to lie eventually along the x-axis 
with footpoints at x — ± 1 . We choose the dimensionless size 
of the rectangular domain tobe !2 x 12 x 8 , a box sufficiently 
large compared to the magnetic structure to avoid boundary 
effects. A variable-size mesh of 127 x 83 x 67 is used. Before 
the emergence of the loop, the initial state consists of a small 
background field as is commonly seen in the solar atmosphere. 
We simulate this background with a simple potential field, spec- 
ified by its normal component at the base to be tanh (x/L) with 
a magnitude of 5% of the expected strength of the loop field. 
The spacetime profiles of the emerging magnetic field and cur- 
rent density are chosen as follows. The two magnetic poles, 
located at the base boundary with opposite polarities, are given 
by B 0 .(x, y, t; R ) = +/ B (r)exp(-[(x - x R ) 2 + r]/a 2 ) for the right 
pole and similarly, with opposite sign, for the left. Here f B {t) 
is the time profile of the emerging field, and (x R (t), 0 ) and 
(x L (t), 0) are the coordinates of the centers of the poles. They 
are initially located at x RM = ± 0.2 at the beginning and move 
to the final locations atx BX = ± 1.0 with the same time profile. 
The Gaussian width of the profile is chosen to be a = 0.5. The 
peak value of B 0;R is shown in Figure 1. 

The current density also consists of two opposite poles with 
their centers coinciding with the magnetic poles at all time. 
The spatial profile is chosen in such a way that the integrated 
net current for each pole is zero. Although this condition is 
not necessary, it is convenient for computational purposes and 
can easily be relaxed to consider a more general case in which 
a reverse-current layer would be self- induced (Van Hoven et 
al. 1995). The profile for the right pole, for example, is J^ix, 
y, t\ R) - 1 - r 2 /c 2 ) exp( - r R /b 2 )B 0 ,(x, y, t, R), where 

fj(t) is the model time profile, r R = (jc - x B ) 2 + y 2 , and 
c~ 2 = a~ 2 + b 2 . The width of the current profile is chosen to 
be b = 0.6a. The parameter a (taken to be negative) is the 
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Fig. 2. — Field lines in the loop region (see text) at times ( top to bottom) t 
= 10, 20. 40. and 360. 


nominal ratio of JIB in a force-free field. (The use of a is 
conventional but does not indicate that the independently 
emerging values of B and J are anywhere near force-free.) The 
peak current density at the center of the right pole is also shown 
in Figure 1 . The profile of the left pole is chosen in a similar 
manner. The purpose of specifying the poles in this way is to 
model the normal components of the flux rope as it surfaces 
dynamically through a quiet photosphere. The value of a, cho- 
sen to be —3.3, is the intended limiting value of the force-free 
parameter, evaluated at the pole center, as the loop relaxes to 
the time-asymptotic, equilibrium, force-free state. In the sim- 
ulation described below, the rise time of the magnetic field is 
set to 100 Alfven transit periods, which is approximately 10-20 
minutes for typical coronal parameters. 

3. RESULTS 

As the subphotospherically generated values of B and J 
emerge, they are far from force-free, as they must be in the 
low-/3 corona. Thus, the fields dynamically evolve, by rising 
and expanding, in moving toward this state. In the following, 
we provide the first description of the details of this oft-imag- 
ined process. 

Figure 2 (top) shows the magnetic field lines at time / = 10. 
The magnetic structure has just emerged from the surface and 
is still close to the photosphere. It lacks the physical appearance 
of a “loop” in the conventional sense as we can expect from 
a low-lying structure. The loop is characterized by the torsion 
parameter a*(r; t) = J • BIB 2 , which generally describes the 
local behavior of the field lines in a sheared field. The contours 
of a*(y, z; t = 10) at the midplane (x = 0) are shown in Figure 
3 (top). The local maximum of | or* | can be considered to be 
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Fig. 3. — (Top) Contours of a * at the midplane (jc = 0) at time / = 10. Only 
the central section, where magnetic flux is emerging, is plotted. The entire 
domain is — 6 < x, y < 6, 0 < z < 8. The base plane limits are a*(min) - 
- 10.1 and a*(max) = 1.5. Dashed lines indicate that the contour values are 
negative. The emerging flux does not yet have the appearance of a loop. Notice 
that the scales of each panel are all different to make the illustration more 
legible. (Second from top) Contours of a* at / = 20. a*(min) = -7.5 and 
a*(max) = 2.5. The flux tube begins to rise from the surface. (Third from 
top) Contours of a* at t — 40. a*(min) = —7.4 and a*(max) = 2.1. The 
loop has levitated from the surface. ( Bottom ) Contours of a* at t — 360. 
of*(min) = - 1.8 and a*(max) = 0.9. The loop has relaxed to its equilibrium 
state. 


the loop “axis,” connecting the two poles on the surface. The 
value of a* near the axis is negative because the local J and 
B have been chosen to be in opposite directions; a*(r; t) in- 
creases away from this axis, then changes sign at some distance 
from the axis, before eventually going to zero. The “apex” of 
the loop axis is at an altitude of only 0.035. Field lines are 
drawn around this axis to depict the nature of the^ magnetic 
structure. The moderate pitch angle of the helical field lines 
clearly shows the effect of current flowing along the loop axis. 

Note that the system is far from being force-free at this early 
time. In fact, the loop is dynamically rising and toroidally 
expanding as driven by the increasing J x B force. Hence, 
a*(r) should not be confused with the conventional a in a 
force-free field, although the former asymptotically approaches 
the latter. At the midplane, the absolute value of a*, which is 
approximately equal to J x /B x locally, is larger than the intended 
a of the final state. However, the loop is not MHD unstable 
because the field lines are relatively short, and the most un- 
stable, long-wavelength modes are not accessible to the loop. 
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The plasma flow at the base has both normal and tangential 
components. The existence of a finite normal component in- 
dicates that there is a net inflow into the box. However, the 
upward flow is confined mostly to a neighborhood directly 
under the apex with a maximum velocity less than 0.02t^. 

Figure 2 {second from top) shows the magnetic field lines 
at t = 20. The magnetic structure has taken on the characteristic 
low-lying loop shape. The corresponding contours of a*(y, z) 
are shown in the second panel of Figure 3. The contours of 
the loop part are not circular because the central part is rising 
into the corona at a speed of 0.023v,, somewhat more quickly 
than either side at 8.3 x 10" The dynamic nature of the 
structure at this instant can be seen by reference to Figure 1 , 
showing that both B (): and J {) , are in rapidly rising phases. Field 
lines are drawn around this a* channel as before. Figure 2 
{third from top) shows the magnetic field lines at t = 40, and 
the third panel in Figure 3 shows the corresponding contours 
of of*. The loop has risen to a substantial altitude so that its 
bottom side has truly detached from the photosphere. The main 
current channel, with negative a*, is surrounded by the return 
current (positive a*). The current density is near its peak value 
in time, and the loop axis is still rising at a speed of 0.022^. 

Finally, we advance the equations until t = 360, at which 
time both J 0: and B 0: have reached their respective asymptotic 
values. The loop has now relaxed to its equilibrium state. The 
field lines and a* are shown in Figures 2 and 3 {bottom panels). 
The altitude of the axis of the loop apex is 1.1, and the radius 
of the main current channel is approximately 0.5. The structure 
of the cross section is now qualitatively similar to Figure 4 of 
Van Hoven et al. (1995). Note that the a contours at the base 
in this earlier model computation are not circular, but the flow 
is, by design. In the present case, the a contours at the base 
are circular by design, while the plasma flow reacts to a mag- 
netic force that is not circular. Hence, these two cases are only 
qualitatively the same even when we try to set a at the center 
of the poles to the same value. 

Because of the finite resistivity, a small electric potential <i>, 
though constant in time, is needed at the base to maintain the 
current. The resulting tangential electric field causes a small 
plasma vortex around each pole as indicated by E x BIB 2 with 
a maximum speed of 4 x 10~V In this state, the magnetic 
loop is only approximately force-free, since the integrated (over 
the entire box) value of \J x B\ is just under 8% of the 
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integrated |7 • B \ . This is due to the relatively small Lundquist 
number we use (10 4 ). The electric potential that keeps the cur- 
rent from decaying due to ohmic dissipation induces an E x 
B drift of the plasma. This finite velocity must be balanced by 
a finite J x B force in equation (2) to maintain a steady state. 
In reality, the Lundquist number in the corona is ~10‘\ and 
the loop is closer to being force-free than the computational 
results indicate. 

4. SUMMARY AND DISCUSSION 

We have described the first demonstration of the dynamic 
emergence of a magnetic loop structure from the solar surface. 
We have taken a phenomenological approach by modeling the 
observed growth of the line-of-sight magnetic field and current 
density on the photosphere, in addition to making a number of 
well-accepted theoretical assumptions. 

As the emerging field and current density at the base reach 
their predetermined values, the magnetic loop gradually settles 
into a nearly force-free equilibrium as the kinetic energy of the 
plasma motions and the integrated magnetic force decay to 
vanishingly small values. The internal structure of the loop, 
including the field geometry and the current distribution, is 
qualitatively similar to the magnetic loop exhibited in our pre- 
vious work (Van Hoven et al. 1995), although the loop is 
formed through an entirely different, observationally moti- 
vated, process and driven by a different mechanism. Finally, 
we want to point out that we have chosen a specific set of 
parameters in this simulation. For example, the time profile of 
the current density is carefully chosen so that the emerging 
magnetic-field lines display visible twists as in observations, 
and yet the loop does not cross the ideal-MHD stability bound- 
ary (Van Hoven 1981). There are many physically acceptable 
regions in the parameter space that we have not explored. In 
addition, there is always the question of possible ideal MHD 
or resistive instability, or flaring, when the loop is overdriven. 
We plan to investigate these issues in future studies. 
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The Initiation of Coronal Mass Ejections by Magnetic Shear 
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Theoretical MHD models have shown that magnetic shear, which can be 
induced in coronal magnetic fields by photospheric flows (including 
differential solar rotation), can lead to the destabilization of large-scale 
coronal magnetic arcades and coronal streamers. When the shear exceeds a 
critical threshold, helmet streamers erupt, with characteristics that are generally 
similar to observations of coronal mass ejections (CMEs). These results 
indicate that magnetic shear can be considered to be a candidate for the 
initiation of CMEs. We discuss the concept of magnetic shear in the corona, 
and we describe its role in the energization of the coronal magnetic field. We 
review some theoretical results on the shearing of axisymmetric coronal 
arcades and streamers, and we present preliminary results on the evolution of 
a three-dimensional model of the solar corona at solar minimum, including 
the eruption of magnetic fields that resemble a CME. 


1. INTRODUCTION 

Although CMEs have now been observed for over two 
decades, we still do not known how they are initiated. The 
principal reason for this lack of understanding is that CMEs 
have been observed on the solar limb, where it is difficult to 
correlate them with changes in the photospheric magnetic 
field. Since the coronal magnetic field cannot be measured 
in general, the magnetic structure of coronal streamers prior 
to, and during eruption, is not known. Two promising 
developments in observations of CME initiation have 
occurred recently. First, a significant correlation has been 
shown to exist between erupting filaments observed on the 
solar disk and emergence of new flux [ Feynman and Martin , 
1995], which implies that emerging flux may initiate 
CMEs, since CMEs frequently (but not always) contain an 
erupting filament. Second, interpretations of Yohkoh SXT 
observations have shown that CMEs may be associated 
with “dimming” of the X-ray corona [e.g., Hudson et al., 
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1996]. These observational techniques can relate CMEs to 
conditions in the photosphere and corona, and may help to 
identify the mechanism of CME initiation. 

In order to develop a complete understanding of CMEs, a 
corresponding theoretical study of CME initiation is 
required. Since a CME expels a large amount of mass (up 
to 10l6 g) into the solar wind and liberates a substantial 
amount of thermal energy in the corona, a successful 
theoretical model must demonstrate how this energy can be 
stored in the corona prior to eruption, in addition to 
showing how this energy can be released impulsively. 

Over the last several years, theoretical models of the 
large-scale corona have evolved considerably. Simple 
idealized models have given way to sophisticated models 
which, in the near future, will allow us to follow the three- 
dimensional evolution and eruption of coronal streamers 
that match measured photospheric magnetic fields. Such 
models are useful for developing intuition about the 
properties of coronal streamers, testing hypotheses about 
CME initiation mechanisms, understanding the signatures 
of interplanetary CME observations, and eventually, being 
able to forecast the trajectory of CMEs, a principal goal of 
the National Space Weather Program. 

In this paper we describe the use of theoretical models to 
study the role of magnetic shear in the initiation of CMEs. 
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Our results indicate that, from a theoretical point of view, 
an increase of magnetic shear may be the mechanism by 
which coronal mass ejections are initiated. The role of 
magnetic shear will undoubtedly be clarified in the future by 
more realistic theoretical calculations and with better 
observations of CMEs. 

2. MAGNETIC SHEAR 

The magnetic field is the principal mechanism by which 
energy can be stored in the lower solar corona. The term 
“magnetic shear” has been used loosely to refer to the 
energized state of the magnetic field in the solar corona. To 
illustrate the definition of magnetic shear, consider the 
idealized situation of a force-free field [e.g., Priest , 1982, 
p. 119], a good approximation to the state of strong 
magnetic fields in active regions, in which magnetic forces 
dominate other forces. In this approximation, the magnetic 
field is the only source of energy. The potential magnetic 
field is the lowest energy state for a given flux distribution 
in the photosphere. Therefore, in order to release energy, 
the magnetic field needs to be “energized” above the 
potential state. The magnetic field can be energized by 
being twisted; the electric currents associated with this twist 
provide a source of free energy that can, in principle, be 
released during an eruption. (Within the force-free model, 
the electric currents associated with this twist would flow 
along the magnetic field.) One way in which the coronal 
magnetic field can be twisted is by (non-uniform) 
photospheric flows. Convective flows in the dense 
photosphere tend to move the footpoints of the magnetic 
field lines that penetrate it through the effect of line tying 
[e.g., Raadu, 1972; Einaudi and Van Hoven, 1981; Priest, 
1982]. We refer to the twist introduced by this effect as 
photospheric shearing. Magnetic twist can also arise in the 
corona from the eruption of twisted (current-carrying) fields 
from below the photosphere [e.g., Leka , 1995; Lites et al., 
1995]. The twist that energizes the magnetic field above 
the potential state has been termed magnetic shear. Note 
that magnetic shear includes twisting by photospheric shear 
flows as well as the twist introduced by eruption of twisted 
magnetic flux from below the photosphere. 

In the context of coronal mass ejections, the situation is 
more complicated than in the idealized case of force-free 
magnetic fields, since additional sources of energy, 
including the kinetic energy in the flowing solar wind, 
gravitational potential energy, and thermal plasma energy, 
become important. Nevertheless, the concept of magnetic 
shear can still be used to express the energization of the 
coronal magnetic field. 

In this paper we discuss magnetic shear introduced by the 
effect of photospheric flows only. The important effect of 
increasing magnetic shear through emerging flux is not 
addressed here. We describe how photospheric shear can 
store energy in the coronal magnetic field, and how it can 


lead to eruption. Such photospheric flows can arise from 
differential rotation and other large-scale flows in the 
photosphere. We therefore refer to the “shearing” of the 
large-scale coronal field by photospheric flows. 

Magnetic field observations that are readily available 
today (e.g., synoptic maps of the normal component of the 
magnetic field in the photosphere deduced from line-of-sight 
magnetograms, from Wilcox Solar Observatory and the 
National Solar Observatory at Kitt Peak) cannot specify the 
magnetic shear in the corona. In order to fully specify the 
state of the corona, and in particular, its level of 
energization, it is necessary to specify the transverse 
component of the magnetic field. Thus, a fundamental 
aspect of the state of solar magnetic field is not provided by 
these magnetograms, as they do not allow configurations 
with different levels of magnetic shear to be distinguished. 
As is well known from studies of magnetic fields in active 
regions, a vector magnetogram is required to uniquely 
specify the magnetic field [Mikic and McClymont, 1994; 
Mikic et al., 1996]. In principle, full-disk vector 
magnetograms can provide information about the transverse 
component of the magnetic field. Whether this can be done 
with sufficient accuracy to determine the relatively weak 
large-scale field is not known at present. In the future, the 
development of this capability will improve our ability to 
model and assess the state of magnetic shear in the solar 
corona. 

Let us illustrate this situation by means of an idealized 
example. Consider the case of an axisymmetric 
equilibrium. The equilibrium magnetic field in the corona 
for a given normal magnetic field distribution in the 
photosphere (e.g., that corresponding to a dipole) can be 
found by solving the steady-state MHD equations with 
finite resistivity. We normally find the steady-state 
solution using a relaxation scheme, by integrating the time- 
dependent MHD equations to steady state with a fixed 
normal component of the magnetic field in the photosphere. 
The resulting solution is a “minimum magnetic shear” 
solution (loosely speaking) for the given boundary 
conditions, since it is found by a relaxation procedure in the 
presence of plasma resistivity. It will be the solution with 
the smallest twist in the magnetic field compatible with the 
specified normal magnetic field. The solution consists of 
the canonical coronal streamer configuration described by 
Pneuman and Kopp [1971], containing a dense closed-field 
region, in which the solar wind flow is arrested, surrounded 
by open magnetic field lines along which the solar wind 
flow streams to supersonic velocities. The magnetic field is 
potential (current-free) nearly everywhere, except in a narrow 
layer surrounding the open-closed field boundary, and along 
a sheet that extends from the tip of the coronal streamer to 
infinity, at which the current is concentrated into a sheet. 
Figure 1 shows an example of a configuration we computed 
with a dipole magnetic flux distribution [Linker and Mikic, 
1995]. The current sheet that borders the open/closed field 
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Figure 1 . (a) Magnetic field lines (contours of constant 

flux y/), and (b) contours of the azimuthal current density (J§) 
in an axisymmetric helmet streamer. A current sheet separates 
fields of opposite polarity and bounds the closed field region; 
away form the current sheet the field is nearly potential. The 
solar wind plasma flows outward along field lines in the open- 
field region, but is essentially stationary inside the closed-field 
region. 


This equilibrium can be energized by twisting the 
magnetic field. In general, photospheric flows will twist 
the magnetic field, introducing magnetic shear. As noted 
above, without detailed observations of the transverse 
magnetic field it is not possible to uniquely determine the 
magnetic shear in the coronal field, so that at present we are 
limited to using phenomenological photospheric shear flow 
profiles to increase the twist of the magnetic field to a level 
that will cause eruption. 

Therefore, the magnetic shear in our “initial” equilibrium 
state is not likely to correspond to that of the corona at any 
time, but is merely a starting point for our computations. 
In particular, this is the reason that a significant amount of 
shear is needed to cause the first eruption in our simulations 
(described below), since we start with a state which has 
minimum magnetic shear. Subsequent eruptions do not 
require as much shear. In principle, once vector magnetic 
field measurements become available, it may be possible to 
construct a coronal magnetic field equilibrium with a level 
of magnetic shear that corresponds to the state of the corona 
corresponding to the particular time of the magnetic field 
observation. At present, lacking more complete 
observational information about the state of the magnetic 
field, our best alternative is to energize the field (by 
applying phenomenological photospheric flows) to a level 
of magnetic shear that causes eruption to occur. 

3. AN MHD MODEL OF THE SOLAR CORONA 


A self-consistent description of the large-scale solar 
corona requires the coupled interaction of magnetic, plasma, 
and solar gravity forces, including the effect of the solar 
wind. In the magnetohydrodynamic (MHD) model, the 
coronal plasma is described by the following equations: 


V xB = — J , 
c 

„ 1 3B 

VxE = - - t , 

c dt 


E + -vxB = 7? J , 
c 1 


( 1 ) 

( 2 ) 

( 3 ) 


boundary is caused by the discontinuity of the magnetic 
field pressure, which in turn is induced by the discontinuity 
in plasma pressure between the open and closed field 
regions. 

In this idealized axisymmetric case, the field lines have an 
arcade-like geometry in the closed-field region, where the 
field is entirely poloidal, with zero azimuthal (B<p) magnetic 
field. In addition, the parallel component of the electric 
current density is everywhere zero. It is this component 
that is associated with the twist in the magnetic field. 
Therefore, this state is not likely to represent the state of 
the corona at any time, since it corresponds to a magnetic 
field that has no twist in it. This field is not likely to erupt 
without additional energization. 


&+V-(pv) = 0, (4) 

p(j^+ v-Vv)=^-JxB-V/>- Vp w 

+ pg + V-(vpVv) , (5) 

V-(pv) = (y- 1)(- pVv + S) , (6) 

where B is the magnetic field intensity, J is the electric 
current density, E is the electric field, v,p, and p are the 
plasma velocity, pressure, and mass density, g is the 
gravitational acceleration, 77 is the plasma resistivity, v is 
the kinematic plasma viscosity, S represents energy source 
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terms, and the wave pressure p w represents the acceleration 
due to Alfven waves. 

The application of this model to the structure and 
dynamics of the solar corona and inner heliosphere is 
discussed by Linker and Mikid [1997, in this volume]. Our 
thermodynamic model is presently being improved by 
including the effects of coronal heating, parallel thermal 
conduction, radiation loss, and acceleration due to Alfven 
waves. We do not consider these improvements here. In 
the present application, we use the poly tropic model 
[ Parker , 1963] with an adiabatic energy equation ( S = 0) 
with 7 =1.05, and no Alfvdn waves (p w = 0). 

We have developed a three-dimensional code to solve 
equations ( 1 )— ( 6 ) in spherical coordinates, as described by 
Mikic and Linker [1994]. This time-dependent model has 
been used to study the evolution of axisymmetric magnetic 
arcades [Mikic and Linker , 1994] and helmet streamers 
[. Linker and Mikic , 1995], as well as the three-dimensional 
structure of the corona [Mikic and Linker , 1996; Linker et 
al. , 1996; Linker and Mikic, 1997]. 

4. DISRUPTION OF AXISYMMETRIC CORONAL 
ARCADES AND STREAMERS 

Our efforts to understand CMEs have focused on first 
distilling the essential physics from the simplest model 
possible (disruption of axisymmetric coronal magnetic 
arcades), and then incorporating these effects into more 
realistic models of the solar corona (including the effect of 
the solar wind, differential rotation, and three-dimensional 
geometry). 

The properties and stability of coronal magnetic fields 
have been studied extensively [see the references in Mikic 
and Linker , 1994]. In order to study the theoretical aspects 
of CME initiation, we started with the simplest model 
possible: we assumed zero beta [i.e., magnetic forces 
dominate plasma forces, so that we can neglect V/? in 
Eq. (5)], a fixed density, we neglected gravity, and we 
modeled two-dimensional (axisymmetric) variation. We 
investigated the dynamical evolution of an initially dipolar 
magnetic field arcade subjected to idealized photospheric 
shearing motions [Mikic and Linker, 1994]. The 
calculations were performed using both the ideal and 
resistive MHD equations. When an arcade is subjected to a 
photospheric shear flow profile, the arcade evolves quasi- 

statically for small amounts of shear. During this phase, 
poloidal magnetic field ( B r , Bq) is converted into azimuthal 
magnetic field (Bq), and the magnetic energy increases with 
increasing magnetic shear as the field becomes more and 
more twisted; energy is therefore stored in the magnetic 
field. However, when a critical shear is exceeded, the field 
expands rapidly and produces a concentration of the electric 
current density. An ideal MHD (i.e., zero resistivity) 
calculation shows that a transition to a partially open 
configuration occurs at the critical shear value. In this state 


a small fraction of the magnetic field lines are closed but the 
majority of field lines (97% of the flux) are open. The open 
field lines of opposite polarity are separated by a tangential 
discontinuity. The magnetic energy of this partially open 
configuration is close to but less than the energy in a fully 
open field [Aly, 1984, 1991; Sturrock , 1991]. 

The transition to a partially open field requires an initially 
smooth magnetic field to evolve into one with 
discontinuities; this process has been described as magnetic 
nonequilibrium [ Parker , 1972, 1979; Priest, 1981; 
Vainshtein and Parker, 1986]. The appearance of a 
discontinuity implies that even a small amount of plasma 
resistivity is important. When we included finite 
resistivity, the discontinuity was resolved into a current 
sheet which was subsequently the site of rapid magnetic 
reconnection, leading to fast flows and the ejection of a 
plasmoid [Mikic and Linker, 1994]. These results suggest 
that CMEs may be initiated by the destabilization of 
magnetic arcades by photospheric shear. The rapid inflation 
of the field at the critical shear value was also found by 
Roumeliotis et al. [1994]. 

Having found the underlying cause of the disruption of 
magnetic configurations, we applied our model to a more 
realistic equilibrium. A comparison with CME 
observations requires the important effect of the solar wind 
to be included. We started with an axisymmetric helmet 
streamer equilibrium corresponding to a dipole magnetic 
field distribution at the solar surface [Pneuman and Kopp, 
1971], which was developed by integrating the time- 
dependent resistive MHD equations to steady state 
[Steinolfson et al., 1982; Linker et al., 1990; Wang et al., 
1993]. The evolution of this field in response to 
photospheric shear flow is qualitatively similar to that of 
the arcade. The closed-field region initially expands slowly 
as the field evolves quasi-statically; when a critical shear is 
reached, the magnetic field lines erupt outward, driving 
plasma into the outer corona. While the underlying reason 
for the disruption is the same in both cases (ideal MHD 
magnetic nonequilibrium), once the streamer begins to rise, 
the plasma within it accelerates into the solar wind, 
stretching and opening the magnetic field lines, and creating 
a current sheet at which the low-lying loops subsequently 
reconnect [Linker and Mikic, 1995]. As the reconnection 
proceeds, the closed-field region grows in size as 
successively higher loops reconnect [Kopp and Pneuman, 
1976], a phenomenon that has been observed in Yohkoh 
soft X-ray images [Hiei et al., 1993; Tsuneta, 1996]. 

These studies of arcades and helmet streamers used 
idealized photospheric flow profiles to induce the magnetic 
shear. One component of photospheric motions that may 
contribute to the energization of large-scale coronal fields is 
differential rotation. We have investigated how differential 
rotation [Snodgrass, 1983] affects axisymmetric coronal 
streamers over many rotations. As in the case of the 
idealized shear profile, the streamer disrupts when a critical 
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Central Meridian: 0 degrees 



Central Meridian: 90 degrees 


Figure 2. An idealized three-dimensional model of the solar minimum corona. Shown on the left are traces of 
magnetic field lines; on the right is the corresponding computed plane-of-sky polarization brightness, (a) The 
view at a central meridian of 0, and (b) the view at a central meridian of 90°. It is apparent that the helmet streamer 
belt is tilted and warped. 


shear is exceeded. With continued differential rotation, the 
streamer disrupts recurrently [Linker et al., 1994], A more 
complete understanding of the role of differential rotation in 
initiating CMEs requires three-dimensional calculations 
with more realistic fields, as described below. 


5. DISRUPTION OF THREE-DIMENSIONAL 
CORONAL STREAMERS 

The axisymmetric results provide, at best, a qualitative 
argument that magnetic shear may be a plausible 
mechanism for the initiation of CMEs. To proceed beyond 


this qualitative agreement with the properties of CMEs, it 
will be necessary to compare the details of the eruption with 
observations (e.g., frequency of eruption, requirements on 
the photospheric shear profile, the signature of the eruption 
in coronagraphs, the properties of the plasmoid in 
interplanetary space). We have already begun this task by 
studying the propagation of an erupted three-dimensional 
plasmoid through interplanetary space to 1 A.U. [Linker and 
Mikic, 1997]. Here we present preliminary results on the 
extension of our studies to the shearing and eruption of 
three-dimensional helmet streamers. 

Our first studies of the evolution and stability of three- 
dimensional magnetic arcades were made using a zero-beta 
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B r in the Photosphere 



Figure 3. The evolution of the radial component of the 
magnetic field in the photosphere, B r , as a function of the 
maximum footpoint displacement, As max . Note that the 
photospheric flow changes the magnetic flux in the solar 
surface. 

model without the solar wind. The results indicated that 
arcades erupt beyond a critical shear, as in the axisym metric 
case. Here we describe the extension of this model to 
include the effect of the solar wind. To generate a 3D 
helmet streamer equilibrium, we first find the potential field 
in the corona that corresponds to the following radial 
magnetic field in the photosphere: 

B r - B 0 {A 0 cos 3 # + A 2 sin 2 # cos20 + A 3 sin 3 # sin30), 

with # 0 = 13.3 Gauss, A 0 = l, A 2 = 0.1, and 
A 3 = 0.025, rotated by -20° about the y-axis. Here (r,0,0) 


are spherical coordinates: # is the co-latitude (# = 0 is the 
North pole; # = 180° is the South pole), and 0 is the 
longitude. (The y-axis is the line # = 90°, 0 = 90° .) This 
is an idealization of the magnetic field at solar minimum, 
with a warped, asymmetric neutral line and relatively 
smooth fields, that resembles the large-scale magnetic field 
for the observational data we have studied, producing a 
tilted, warped heliospheric current sheet. 

The potential field corresponding to this normal field at 
the photosphere was computed as an initial state, and the 
resistive MHD equations were integrated to steady state to 
compute a 3D coronal streamer equilibrium. Figure 2 
shows the magnetic field lines in the equilibrium, along 
with the computed polarization brightness, at two choices 
of central meridian longitude. It is readily apparent that this 
is a three-dimensional configuration. 

The evolution of this equilibrium was then followed in 
response to an applied photospheric shear flow. We used an 
idealized shear profile similar to that used previously for 
axisymmetric arcades [Mikic and Linker , 1994], with 
v = v 0 (#)$, and a width A# m = 30°. (Note that this 
profile does not correspond to differential rotation; the 
evolution of this field for a differential rotation profile is 
discussed by Linker and Miki 6 [1997].) In this three- 
dimensional equilibrium, the normal component of the 
magnetic field in the photosphere changes as a result of the 
advection of magnetic flux. Figure 3 shows the evolution 
of B r in the photosphere. We parameterize the amount of 
shear introduced at any particular time in terms of the 
maximum displacement of a field line footpoint from its 
initial position, Aj max , as in [Mikic and Linker, 1994]. 

In order to minimize the computational time required to 
perform this 3D numerical simulation, we used a shear flow 
velocity that is ~ 10 times larger than flows that are 
typically observed in the photosphere. The maximum shear 
flow velocity was 4.8 km/s, compared to typical 
photospheric flow velocities of 0.5-1 km/s; the 
photospheric flow velocity associated with solar rotation is 
2 km/s. Our previous results [ Mikic and Linker, 1994] 
imply that magnetic fields evolve quasi-statically prior to 
eruption. Our enhanced shear flow velocity is therefore 
expected to shorten the time required to reach a level of 
shear that leads to an eruption, and is not expected to affect 
the nature of the eruption significantly (as long as the shear 
flow velocity is small compared to the Alfven speed). It is 
possible to use a physical value of the shear flow velocity 
in our simulations (at greater computational expense), and 
this will be done in future refinements of this calculation. 

Figure 4 shows the evolution of selected magnetic field 
lines as a function of the footpoint displacement. Note that 
the field lines rise slowly as they become twisted. When 
the shear reaches a critical value (As max ~ 1.6/? 0 , where 
R 0 is the solar radius), the field lines begin to expand 
outward rapidly. This behavior is qualitatively similar to 
that observed in the axisymmetric case, although the field 
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Field Line Evolution 






Figure 4. The evolution of selected magnetic field lines in a three-dimensional configuration as a function of 
the maximum footpoint displacement, A,y max . The field lines initially evolve quasi-statically as they become 
twisted by the applied photospheric shear flow. However, when the twist approaches A.s max ~ 1.6/? 0 , there is a 
rapid upward motion of the field lines. This eruption is characteristic of a coronal mass ejection. 


line geometry is considerably more complicated. The 
erupting magnetic field lines have a long wavelength 
structure, indicating that the eruption is global, with a large 
longitudinal extent. This feature of the eruption agrees with 
recent LASCO coronagraph images of CMEs, in which 


CME ejecta appear to leave the East and West solar limbs 
simultaneously. 

The magnetic energy of the configuration decreases 
following the eruption, as magnetic reconnection takes 
place. The nature of the reconnection is complicated, and is 
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not easily diagnosed. (We deduce that reconnection is 
taking place because the magnetic footpoint connectivity in 
the photosphere changes rapidly after the eruption.) The 
magnetic flux that is ejected has been followed as it 
propagates into interplanetary space in a related simulation 
[Linker and Mikic , 1997]. 


6. DISCUSSION 

We have studied the evolution of 2D and 3D helmet 
streamer configurations in the presence of photospheric 
shearing flows. The streamers initially evolve quasi- 
statically, causing the closed field region to grow in size. 
When a critical shear is reached, the configurations erupt, 
sending plasma and magnetic fields into the outer corona. 
This eruption appears to be an ideal process, related to the 
tendency of the field to open when the energy of the field 
approaches the open field energy. Subsequent to this 
outward eruption of the field lines, reconnection of the 
magnetic field occurs. The nature of this reconnection 
process in three dimensions is quite complicated and is 
presently under investigation. 

Our results indicate that magnetic shear may initiate 
coronal mass ejections. At present, this evidence must be 
regarded as qualitative, since it is based on the study of 
idealized magnetic field configurations. In future work it 
will be necessary to compare the predictions from our 
theoretical model with available observational data, 
especially the properties of plasmoids and “flux ropes” 
observed in the interplanetary medium, to obtain a more 
quantitative evaluation of the role of magnetic shear in 
CME initiation. 
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Extending Coronal Models to Earth Orbit 
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Solar wind conditions at Earth play a primary role in the initiation of ge- 
omagnetic activity. The forecasting of solar wind conditions at Earth based 
on remote observations of the Sun is thus a key element of space weather 
prediction. We describe how observations of the photospheric magnetic field 
can be incorporated into three-dimensional MHD computations of the solar 
corona and inner heliosphere. We show that the resulting solutions com- 
pare favorably with observations and that this same capability can be used 
to model the initiation of coronal mass ejections and their propagation out 
into the heliosphere. These encouraging results suggest that an operational 
computational solar wind model can eventually be developed, suitable for 
forecasting solar wind properties at Earth. 


1. INTRODUCTION 

Coronal mass ejections (CMEs) are exceedingly com- 
plex phenomena. From their initiation on the Sun to 
their propagation through the heliosphere, CMEs span 
a large range of both distance and physical parameter 
space. Understanding how CMEs, typically observed 
as loop-like structures in white-light coronagraphs [e.g., 
Hundhausen, 1993], are initiated and how they ultimate- 
ly manifest themselves in interplanetary space is a fun- 
damental challenge for solar and heliospheric science. 
Apart from the intellectual attraction of such a chal- 
lenge, solution of this problem has significant practical 
applications. It is well known that solar wind condi- 
tions upstream of the Earth’s magnetosphere play an 
important role in geomagnetic activity, and that CMEs 
in particular are associated with the largest geomagnetic 
storms [Gosling, 1993]. Geomagnetic storms can cause 
disruption of satellite operations, communications, nav- 
igation, and electric power distribution grids, and cre- 
ate a hazardous environment for astronauts engaged in 
extra- vehicular activities. The prediction of such “space 
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weather” phenomena has thus become recognized as an 
important problem for the space science community, 
as evidenced by the National Space Weather Program 
Strategic Plan [ Wright et al, 1995]. 

Remote observations of the magnetic and plasma en- 
vironment of the Sun have been made routinely for some 
time. With the SOHO spacecraft now operational, the 
amount and quality of such measurements, including the 
detection of CMEs, has greatly increased. In the context 
of space-weather forecasting, using remote solar obser- 
vations to accurately predict the characteristics of the 
solar wind at Earth orbit, especially the arrival time and 
properties of CMEs, is one of the primary services to be 
provided by solar and heliospheric science. 

To predict effects at Earth from events occurring on 
the Sun, solar observations must be incorporated into 
a physical model. The magnetohydrodynamic (MHD) 
fluid description is an appropriate starting point for 
modeling the solar wind. Even with multi-fluid and ki- 
netic effects neglected, consideration of the important 
physical processes in multi-dimensional geometry ren- 
ders the MHD equations intractable to an analytic ap- 
proach, and a computational solar wind model is neces- 
sary if we are to forecast solar wind conditions at Earth 
orbit. This is not surprising, as computational mod- 
els of the atmosphere have long played an important 
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role in terrestrial weather prediction [Houghton, 1977]. 
Our computational solar wind model must fulfill two 
requirements: (1) computation of the “background” so- 
lar corona and solar wind, and (2) calculation of the 
initiation and propagation of CMEs. Requirement (1) 
arises because the geoeffectiveness of CMEs is related to 
the structure of the interplanetary magnetic field (IMF) 
[Gosling et ai, 1990]; [Crooker et al., 1992]. Also, apart 
from the effect of CMEs, the background solar wind 
plays an important role in geomagnetic activity. It has 
long been known that 27 day recurrences in geomag- 
netic activity are directly linked to the solar rotation 
period [e.g., Maunder, 1905]. Fast streams in the solar 
wind, which originate in coronal holes [e.g., Altschuler 
et ai, 1972], are generally believed to be the cause of 
recurrent geomagnetic storms [Hundhausen, 1979; Zir- 
in, 1988; Foukal , 1990]. Recently the role of the streamer 
belt and corotating interaction regions in producing re- 
current geomagnetic activity has also been recognized 
[ Crooker and Oliver, 1994]. 

At the outset, we must recognize the formidable na- 
ture of the goal we have described and how distant we 
are from achieving it. While many obstacles contribute 
to the difficulty of solar wind forecasting, perhaps the 
most obvious problem is our lack of understanding of 
how CMEs are initiated. Nevertheless, to make progress 
on this task, we must outline a strategy for model de- 
velopment. With this purpose in mind, we demonstrate 
in this paper how MHD models of the solar wind can 
supply the two key requirements for solar wind forecast- 
ing. Section 2 briefly discusses the methodology of our 
computations. In section 3 we describe realistic compu- 
tations of the solar corona, and we show that our results 
compare favorably with coronal and heliospheric obser- 
vations. In sections 4 and 5 we demonstrate a model 
computation of the solar corona and inner heliosphere 
at solar minimum, and we show how one candidate pro- 
cess for CME initiation (shearing of the magnetic foot- 
points by differential rotation) can be modeled and the 
heliospheric effects of the resulting disturbance studied. 
Section 6 summarizes our present capabilities and indi- 
cates future directions. 

2. METHODOLOGY 

To compute MHD solutions for the large-scale corona, 
we solve the following form of the equations in spherical 
coordinates: 

4 TT 

V x B = — J (I) 
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1 dB 

c dt 


-V x E 


E + 


v x B 
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= 77 J 


! + V.(,v) = 0 
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= - J x B — Vp — Vp w 
c 

+ />g + V • (i/pVv) (5) 

^ + V.(pv) = ( 7 — 1 )(-pVv + 5) (6) 

where B is the magnetic field intensity, J is the electric 
current density, E is the electric field, v, p, and p are the 
plasma velocity, mass density, and pressure. The gravi- 
tational acceleration is g, 7 is the ratio of specific heats, 
r) is the resistivity, v is the viscosity, S represents energy 
source terms, and the wave pressure p w represents the 
acceleration due to Alfven waves [ Jacques , 1977; Hoi - 
weg, 1978]. 

The term S in equation (6) includes the effects of 
coronal heating, thermal conduction parallel to B, ra- 
diative losses, and Alfven wave dissipation (viscous and 
resistive dissipation can also be included). A simplified 
model of the corona, known as the “polytropic model” , 
is obtained when an adiabatic energy equation with a 
reduced 7 is used [Parker, 1963]. This is a crude way of 
modeling the complicated thermodynamics in the coro- 
na with a simple energy equation. This choice results 
by setting S — 0 in Eq. (6) and pw — 0 in Eq. (5). With 
this model, values of 7 close to 1 (7 = 1.05 for the re- 
sults shown in this paper) are necessary to produce den- 
sity and temperature profiles that are similar to coronal 
observations; this indicates that the terms included in 
S are in fact important for describing the energy bal- 
ance of the corona. In this paper we describe computa- 
tions using the polytropic model. Computations using 
the full equations (1-6) have been performed [Mikic et 
al., 1996ab] and will be described in a future paper. 

Mikic and Linker [1994] describe the method used to 
solve equations (1-6). To compute coronal and helio- 
spheric solutions, the equations must be supplemented 
with appropriate boundary and initial conditions. In 
spherical geometry, two boundaries appear in the simu- 
lation: the physical inner radial boundary at r = R s (the 
solar radius) and an artificial outer radial boundary at r 
= Ri, which we typically place in the range 20-215 R s . 


dv 

5T + V 


Vv 
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Specification of the appropriate boundary conditions is 
facilitated by examining the characteristic form of equa- 
tions (1-6) [Courani and Friedrichs , 1948]. (See Hu and 
Wu, [1984] for an example using the MHD equations.) 
Characteristics traveling into the domain require that 
physical information be provided. At the inner bound- 
ary, four characteristics point into the domain, and re- 
quire that four quantities be specified. We specify the 
distributions of p and p at r = R s . When no surface 
motions are included in the calculation, we also specify 
that Eg and — 0 at r = R s ; this condition fixes the 
radial magnetic field at the inner boundary (5 r p) at its 
initial value and is equivalent to setting vj_ (velocity 
perpendicular to B) = 0 there. When surface motions 
such as the solar rotation are included, Eg and E^ are 
specified to be consistent with this motion (from equa- 
tion 3); note that in this case the distribution of B t q can 
be modified by the surface motions. The outer bound- 
ary is typically placed well beyond the critical points, so 
all characteristics are outgoing and no physical bound- 
ary conditions are required. To advance the solution at 
the outer boundary, we use the characteristic equations 
to compute v. The plasma (3 (= Srrp/B^ , the ratio of 
plasma pressure to magnetic pressure) is typically 1 or 
greater at the outer boundary, and we find that charac- 
teristics based on the gas equations are sufficient. The 
staggering of the mesh then allows all other quantities to 
be computed in the same manner as the interior points 
[Mikic and Linker , 1994]. Characteristic equations are 
also used to compute vy (velocity parallel to B) at the 
lower boundary. 

For the initial condition, a potential magnetic field 
consistent with the specified distribution of B r at the 
lower boundary, and a wind solution [ Parker , 1963] con- 
sistent with the specified p and p are used. Equations 
(1-6) are then integrated forward in time until a steady 
state is reached. The computations are performed on 
a mesh that is nonuniform in the r and 6 directions: 
A r « .01 R$ near the inner boundary and A r & 10 
R 5 at r = 215 R 5 ; A# varied between .03 and .06 ra- 
dians. The longitudinal (<f>) coordinate is treated using 
a pseudospectral method (this requires a uniform dis- 
tribution of points in <j>). Our higher resolution cases 
used 101 x 101 x 64 (r points; cases extending out 
to approximately 1 A.U. (1 astronomical unit = 1.49 x 
10 6 km = 214 solar radii) used 111 x 51 x 32 points. 

Previous coronal and solar wind solutions of (1-6) 
have typically been performed with idealized magnetic 
fields [Endler, 1971; Pneuman and Kopp, 1971; Steinolf- 


son et al., 1982; Washimi et al ., 1987; Linker ei ai, 1990; 
Wang ei ai, 1993; Linker and Mikic, 1995], or with an 
inner boundary beyond the critical points [Smith and 
Dryer 1990; Detman ei ai, 1991; Pizzo, 1991; Odstr- 
cil, 1994]. To perform a realistic 3-D MHD computa- 
tion of the corona that can be compared with specific 
observations, it is necessary to incorporate solar obser- 
vations into the boundary conditions [Usmanov, 1993; 
Mikic and Linker, 1996; Linker et ai, 1996]. One of 
the most readily available observational data sets is the 
magnetic field at the photosphere. This is also the most 
important observation to address for coronal and helio- 
spheric modeling. We have used Wilcox Solar Observa- 
tory synoptic maps (collected during a solar rotation by 
daily measurements of the line-of-sight magnetic field at 
central meridian) to specify the radial magnetic field at 
the photosphere (in the manner described by Wang and 
Sheeley [1992]). 

3. COMPARISONS WITH CORONAL 
AND HELIOSPHERIC DATA 

Solutions obtained in the manner described in section 
2 can in principle provide a 3D description of the corona 
and inner heliosphere, including the detailed distribu- 
tion of magnetic fields, currents, plasma density, and 
temperature. However, the validity of this approach can 
only be verified through comparison with observations. 
As a test of our coronal modeling capability, we used 
our computations to predict the large-scale structure of 
the solar corona during the October 24, 1995 eclipse 
(occurring during Carrington rotation (CR) 1901), visi- 
ble in a number of locations in the eastern hemisphere. 
We carried out a simulation using photospheric mag- 
netic field data from the previous rotation (CR1900; 
September 2 - September 29, 1995) on October 5, 1995, 
and put the results on the World Wide Web (more de- 
tailed comparisons and new results can be viewed at 
http:/ /iris023.saic.com:8000/corona/modeling.html). 

We also presented the results at the Sacramento Peak 
workshop on October 18, 1995 [Linker ei ai, 1996]. Fig- 
ure 1 (leftmost frame) shows the magnetic field lines 
from our calculation. The view angle corresponds to the 
approximate time of the eclipse. The solution shows the 
formation of helmet streamers; these are regions with 
closed magnetic fields that trap coronal plasma flowing 
out of the Sun. Along open magnetic field lines, the 
solar wind streams freely, reaching supersonic speeds. 

To directly compare our results with observations, we 
develop images of the polarization brightness (pB\ pro- 
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Figure 1. A prediction of the structure of the solar corona during the October 24, 1995 solar eclipse. The MHD 
simulation was carried out on October 5, 1995, using Wilcox synoptic magnetic data for the previous rotation. From left 
to right, the frames show: Field lines, polarization brightness computed from the simulation, and an eclipse photograph 
taken by F. Diego (UCL) in white light with F=910 mm and a two-second exposure time. 


portional to the line-of-sight integral of the product of 
the electron density and a scattering function that varies 
along the line of sight). This quantity is frequently ob- 
served with coronagraphs. Using the plasma density 
from our coronal model, we can compute pB to simulate 
an eclipse or coronagraph image and compare it with the 
actual data. Radially graded filters are applied to eclipse 
images to compensate for the rapid fall-off of coronal 
density with radial distance; we detrend our computed 
pB in a similar manner. The polarization brightness 
of the corona predicted by our simulation, as it would 
be seen on October 24, 1995, at 05:00UT is shown in 
Figure 1 (middle frame), along with an image of the 
eclipse taken by F. Diego of University College, London 
(rightmost frame). The helmet streamers and open field 
regions predicted by the computation agree reasonably 
well with the eclipse observations. We have performed a 
similar comparison for the November 3, 1994 eclipse and 
CR1888 [Mikic and Linker , 1996; Linker ei al. , 1996]. 
These computations support the long-held belief that 
the magnetic field distribution on the Sun controls the 
position and shape of the streamer belt. 

We have also compared the results of our calcula- 
tions with interplanetary observations. As a first test, 
we computed an MHD model of the solar corona for 
CR1869 (May-June 1993). This rotation was of par- 
ticular interest for Ulysses observations, as the Ulysses 
spacecraft ceased to observe sector-boundary crossings 
during that time period [Smith ei al. , 1993]. Figure 
2 shows a comparison of the heliospheric current sheet 
(HCS) predicted by our MHD computation with that 
of the source-surface model [e.g., Schaiien ei al ., 1969; 
Altschuler and Newkirk , 1969; Hoeksema, 1991; Wang 


and Sheeley 1988, 1992], a frequently used tool for ap- 
proximating heliospheric structure. Ulysses’ latitude 
position for this time period (near 30° S latitude) is also 
shown. The source-surface model predicts crossings for 
this time period, whereas the MHD simulation correctly 
predicts no HCS crossings. 

During February-April of 1995 (before and after the 
spacecraft approached perihelion), the Ulysses space- 
craft sampled a wide range of heliographic latitude in 
a short period of time. Figure 3 shows the HCS pre- 
dicted by our MHD computation for CR1892, the start 
of this fast latitude scan. Also shown is the Ulysses tra- 
jectory projected in solar latitude and Carrington longi- 
tude (back at the Sun) and published Ulysses HCS cross- 
ings indicated by crosses [Smith et al., 1995]. The differ- 
ent line styles on the trajectory plots indicate the Car- 
rington rotation at that time. During CR1892 (the time 
period for which the calculation is most valid), the two 
Ulysses crossings occur almost exactly where predicted 
by the MHD computation. Later in time (CR1893 and 
CR1894), the overall shape of the MHD HCS agrees 
well with Smith et al.’s empirically derived HCS, but 
the Ulysses crossings occur above the MHD HCS. The 
reason for this can be seen in Figure 4, which shows the 
predicted source-surface model for the 3 rotations. The 
source-surface model suggests that the solar magnetic 
field is changing during this time period, as evidenced 
by the changing HCS. Therefore, MHD computations of 
CR1893 and CR1894 are required for a complete com- 
parison; this work is presently underway. 

While the favorable comparisons between our com- 
putational results and coronal and heliospheric obser- 
vations are encouraging, it should be noted that there 
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Figure 2. A comparison of the heliospheric current sheet 
predicted by the source-surface model and an MHD calcu- 
lation for Carrington rotation (CR) 1869 (MaylO-June 6, 
1993). The Ulysses spacecraft, which did not observe cur- 
rent sheet crossings during this rotation, was situated at 30° 
latitude. The circles indicate the crossings predicted by the 
source-surface model. 



Heliographic Longitude 


Figure 3. The heliospheric current sheet (HCS) predicted 
by the MHD model for CR1892, with the Ulysses trajectory 
for the fast-latitude scan superimposed. HCS crossings iden- 
tified by Smith et al. [1995] are indicated by black crosses. 
The times of the different rotations (CR1892, CR1893, and 
CR1894) are coded by line style on the trajectory plot. 

are also some differences between our simulations and 
observations. Fine-scale details of the corona do not 
appear in our computations. Higher resolution magne- 
tograms (such as those from the National Solar Obser- 
vatory at Kitt Peak or the SOI/MDI instrument aboard 
SOHO), coupled with higher resolution computations, 
may help to capture some of these fine-scale features. 
Streamers in eclipse images typically show a stronger 
nonradial tendency than in our the computations. This 
may be related to the poor estimation of polar fields 
in the Wilcox data, due to projection effects, and may 
also be improved by better magnetograms. Most impor- 
tant, our computations (using a poly tropic model) fail 
to reproduce the fast (800km/s) solar wind observed by 
Ulysses at high latitude. Improvement of this aspect of 
the calculation requires consideration of the momentum 
and energy source terms discussed in section 2. Our pre- 
liminary ID and 2D calculations including these terms 


Figure 4. Variation of the heliospheric current sheet pre- 
dicted by the source-surface model for the rotations occuring 
during the fast-latitude scan. The extent of the HCS varies 
during this time period. 


show promising results [ Mikic et al ., 1996ab], and fur- 
ther investigation of these solutions is ongoing. 

4. A MODEL OF THE CORONA AND 
INNER HELIOSPHERE AT SOLAR MINIMUM 

To demonstrate how CMEs may be initiated in the 
corona and propagate out into the heliosphere, we devel- 
oped a model configuration of the solar corona (and its 
extension to 1 A.U.) at solar minimum. Guided by our 
experience with the Wilcox photospheric magnetic field 
data and our comparisons with eclipse images, we spec- 
ified an initial magnetic flux distribution of the form: 

B r = Aqcos^^ + Ajsin^0cos20 + A 2 sin^ $sin3(j) (7) 

with Aq = 13.3 G (Gauss), Ai = 1.3 G, A 2 = 0.33 G, 
and the distribution rotated by -20° around the y axis. 
We then computed an equilibrium configuration by in- 
tegrating the MHD equations to steady state, as de- 
scribed in section 2, with the additional constraint that 
the Sun’s rigid rotation rate was imposed (corresponding 
to a sidereal period of 26 days). The resulting config- 
uration is shown in Figure 5. The magnetic field lines 
and polarization brightness near the Sun (Fig. 5a and 
5b) show a configuration similar to that often seen at 
solar minimum. As we move farther from the Sun, the 
magnetic field lines show the expected spiral behavior 
(Figure 5d shows the field lines out to 1 A.U.). With 
this configuration, we can investigate how different pro- 
cesses can affect coronal evolution, how CMEs might 
be initiated, and the heliospheric consequences of these 
events. 
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Figure 5. An MHD simulation of the solar corona and inner heliosphere for a solar-minimum type configuration. The 
computation is performed in the inertial frame, so the magnetic flux distribution on the Sun rotates rigidly, (a) Field lines 
viewed close to the Sun, showing a helmet streamer configuration, (b) polarization brightness from the same view as (a), 
(c) Field lines from the same view angle as (a) and (b), but farther from the Sun. (d) Field lines from 1 A.U. above the 
Sun’s North pole. The spiral structure is apparent. Field lines that appear “shorter” actually are receding from or 
approaching the viewpoint. 


5. MODELING CME INITIATION 
AND PROPAGATION 

While the exact cause of CMEs is unknown, it is gen- 
erally agreed that the amount of energy in the coronal 
magnetic field (above the energy of the corresponding 
potential field) is an important factor in determining if 
and when the coronal magnetic field erupts. Unfortu- 
nately, line-of-sight magnetograms do not contain infor- 
mation about the amount of parallel current that is flow- 
ing in the coronal magnetic field, so we do not know the 
magnetic energy state of the corona from these measure- 
ments. As is well known from studies of active regions, a 
vector magnetogram is required to uniquely specify the 
field. Full-disk vector magnetograms, if available with 


sufficient accuracy, could allow us to compute coronal 
configurations that match the measured twist (or shear) 
in the photospheric magnetic field and determine the 
magnetic energy state of the corona. Lacking this infor- 
mation, the coronal equilibria we compute (such as those 
shown in Figure 1 and 5) correspond to helmet stream- 
ers with a minimum amount of twist in the magnetic 
field and probably have an unrealistically low magnetic 
energy. 

One process that might initiate CMEs is shearing of 
the magnetic field footpoints by photospheric motions 
[Mikic and Linker , 1994; Linker ei al., 1994; Romelio- 
tis , 1994; Linker and Mikic , 1995]. Mikic and Linker 
[this volume] describe how an idealized photospheric 
shear profile causes the coronal configuration of Figure 
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5 to evolve, resulting in eruption of a portion of the 
closed field lines. Because the calculation is starting 
with a magnetic field that is in a much lower energy 
state than the actual corona, an artificially long shear- 
ing time is necessary to initially energize the field. As we 
have shown previously [Linker et ai 1994; Linker and 
Mikic , 1995], after a sheared helmet streamer erupts, 
re-formation of the helmet streamer by magnetic recon- 
nection releases only a portion of the available magnetic 
energy and does not return the configuration to the pre- 
sheared state. Subsequent eruptions thus require less 
shear. 

Here we discuss a computation where a differential 
rotation profile, rather than the idealized shear profile 
described by Mikic and Linker [this volume], was intro- 
duced. To accomplish the initial energization of the field 
more rapidly (and reduce computing time), the rotation 
rate of the Sun was increased by a factor of 10 when we 
introduced differential rotation. The resulting disrup- 
tion of the helmet streamer configuration was similar to 
that described by Mikic and Linker [this volume]. The 
introduction of shear initially causes the magnetic field 
to expand slowly. When a critical shear is reached, the 
magnetic field erupts rapidly outward, followed by re- 
connection of magnetic field lines. In this respect our 
3D results are similar to previous 2D results [Linker et 
at., 1994; Linker and Mikic , 1995], but the evolution 
of the magnetic field in 3D, particularly the magnetic 
reconnection, is more complicated. As an example of 
how we can investigate the heliospheric consequences of 
such eruptions, Figure 6 shows the magnetic field line 
evolution out to 0.5 A.U. (a portion of the total sim- 
ulation domain). The black field lines extend out to 1 
A.U. in the equilibrium shown in Figure 5; note that the 
artificially increased rotation rate of the Sun increases 
the spiral angle for these field lines. A magnetic erup- 
tion at t = 44 hours results in a portion of the helmet 
streamer of Figure 5 expanding rapidly outward into the 
heliosphere (the gray field lines). A subsequent eruption 
beginning at t = 76 hours results in more magnetic flux 
being carried out into the heliosphere. 

In the simpler geometry of the previously mentioned 
2D studies, a completely detached plasmoid (a torus sur- 
rounding the Sun) propagated outward. In the more 
complicated 3D case shown here, no magnetic field lines 
that are completely detached from the Sun are apparent, 
but their presence has not yet been ruled out. These cal- 
culations represent our first efforts to investigate coronal 
mass ejections in three dimensions and should only be 







Figure 6. Evolution of magnetic field lines viewed from 
above the Sun’s north pole at 0.5 A.U., after the solar- 
minimum configuration of Figure 5 is subjected to an en- 
hanced differential rotation profile. Black field lines were 
open in the initial configuration; gray field lines were ini- 
tially part of the closed-field helmet streamers. Magnetic 
eruptions at 44 hours and again at 76 hours cause magnetic 
flux that previously closed near the Sun to be carried out 
into the heliosphere. 

regarded as a first step. We plan in further studies to 
examine the details of the magnetic topology both in 
the corona and far from the Sun, as well as studying 
other possible initiation mechanisms. 

6. SUMMARY 

An important element of predicting geomagnetic ac- 
tivity is forecasting solar wind conditions at Earth or- 
bit. Essential to this effort is a computational model 
of the solar wind capable of describing (1) the struc- 
ture of the solar corona and inner heliosphere, and (2) 
the initiation and propagation of coronal mass ejections. 
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We have described an MHD model that can in princi- 
ple provide these capabilities. The favorable compar- 
isons of our model with coronal and heliospheric data 
indicate that the first of these goals is the most feasi- 
ble, although a more sophisticated treatment of energy 
transport in the solar wind is necessary to accurately 
compute solar wind velocities. The second goal is the 
more difficult task, since our present understanding of 
CMEs and their initiation is limited. Simulations like 
those we have described are an attempt to understand 
the basic phenomena of CMEs. With the SOHO mis- 
sion now operational and with continued Yohkoh and 
ground-based observations, the next few years should 
see a rapid growth in our knowledge of CMEs. The con- 
fluence of improved observations and more sophisticated 
theory and modeling may lead to an improved under- 
standing of CMEs and, eventually, to the capability to 
forecast CME effects at the Earth. 
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